#Crevasses on moons 
#calculates fracture properties with variable h
#robert.law@uib.no 2024

import os
import sys
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.integrate import quad

os.chdir(os.path.dirname(sys.argv[0]))
dir = os.path.dirname(os.path.realpath(__file__))

#equations from https://www.sciencedirect.com/science/article/pii/S0165232X97000220

#variables
H = 30 * 1000 #(m) depth of ice below fracture
h = 0 #(m) depth of ice above fracture
rhoi = 917 #(kg m-3) density of ice
rhos = 917 #(kg m-3) surface density
rhow = 1000 #(kg m-3) density of water
CC = 0.02 #(m-1) constant for density calculations
g = 1.31 #(m s-2) gravity 
Rxx = 0.2*1e6 #(Pa) background tensile stress
d = np.arange(0.5, 400, 1) #(m) crevasse depth
a = 500

nice_fig = True #make the nice figure
save_nice_fig = False#save the nice figure
Ks = False #plot ks plot

#h_range = np.arange(0, 1001, 50)
h_range = np.array([0, 50, 100, 200, 500, 1000])
#a_range = np.arange(0, 1001, 500)

#initialise output
Kic = np.zeros((np.size(d), np.size(h_range)))
Ki1 = np.zeros_like(Kic)
Ki2 = np.zeros_like(Kic)
Ki3 = np.zeros_like(Kic)
flambda_out = np.zeros_like(Kic)



#define G 
def G(b, d, H):  
    
    #Eq. 12
    gamm_a = b/d
    lambd_a = (d+h+a)/H

    A = ( 3.52*(1-gamm_a) )/( (1 - lambd_a)**(3/2) )
    B = ( 4.35 - 5.28*(gamm_a) )/( (1 - lambd_a)**(1/2) )
    C = ( 1.30 - 0.30*( (gamm_a)**(3/2 ) ) )/( (1 - ( (gamm_a)**2 ))**(1/2) )
    D = 0.83 - 1.76*(gamm_a)
    E = 1 - (1 - gamm_a)*(lambd_a)

    return A - B + (C + D)*E


#for loop 
for j in range(len(h_range)):
  print('a is: ' + str(h_range[j]) + ' m')
  for i in range(len(d)):

    #print('Depth is: ' + str(d[i]) + ' m')

    #fix functions to values of d and H required
    #int2_fixed = lambda b: (b + ((rhoi -rhos)/(rhoi*CC))*(1 - np.exp(-CC*b)) )*G(b, d[i], H)
    int2_fixed = lambda b: (-b*rhoi - h_range[j]*rhoi - a*rhow)*G(b, d[i], H)
    int3_fixed = lambda b: (b + a)*G(b, d[i], H)

    #calculate integrals within bounds
    int2, err = quad(int2_fixed, 0, d[i])
    #if a > 0:
    #  int3, err = quad(int3_fixed, a, d[i])
    #else:
    int3, err = quad(int3_fixed, 0, d[i]) #to account for body of overlying water

    #calculate F(lambda) Eq. 6
    lambd_a = (d[i] + h_range[j] + a)/H
    flambda = 1.12 - 0.23*((lambd_a)) + 10.55*((lambd_a)**2) - 21.72*((lambd_a)**3) + 30.39*((lambd_a)**4)
    flambda_out[i, j] = flambda

    #calculate stress intensity (Eqs. 5, 14, 16)
    Ki1[i, j] = flambda*Rxx*( ( np.pi * d[i] )**(1/2) )
    Ki2[i, j] = ( ( 2*g )/( ( np.pi * d[i])**(1/2) ) ) * int2 #- (rhoi*g*h)/ (( np.pi * d[i] )**(1/2))
    Ki3[i, j] = ( ( 2*rhow*g )/( ( np.pi * d[i])**(1/2) ) ) * int3
    #if d[i] < a:
    # Ki3 = 0 #set Ki3 to 0 if crevasse depth is less than prescribed water level 

    Kic[i, j] = Ki1[i, j] + np.nan_to_num(Ki2[i, j]) + np.nan_to_num(Ki3[i, j])

if nice_fig:
   
  mpl.rcParams['font.size'] = 8
  plt.figure(figsize=(80/25.4, 80/25.4))

  for j in range(len(h_range)):
    color = plt.cm.Greys_r(j  / (len(h_range) + 4))
    plt.plot(d, Kic[:,j]/1e6, color=color, lw = 1)

  plt.plot(d, Kic[:,0]/1e6, lw=2, color='black')
  
  plt.xlabel('Fracture depth ($m$)')
  plt.ylabel('Stress intensity factor ($MPa^{1/2}$)')
  plt.axhline(y=0.4, color=(0.76, 0.83, 0.89, 1), linestyle='--', lw=2)
  #plt.axhline(y=0, color='black', linestyle='-', lw=0.5)
  plt.xlim(0, 200)
  plt.ylim(-25, 10)
  plt.tight_layout()

  if save_nice_fig:
    plt.savefig(f'Fracture_depth_Rxx_{round(Rxx*1e-6,2)}.png', dpi = 400)

  plt.show()

  if Ks:
    #plot K1, K2, K3

    plt.figure(figsize=(80/25.4, 80/25.4))
    
    plt.plot(d, Ki1/1e6, lw = 1, color=(0.6667, 0.0235, 0.9608, 1))
    plt.plot(d, Ki2/1e6, lw = 1, linestyle = '--', color=(0.3373, 0.4039, 0.8706, 1))
    plt.plot(d, Ki3/1e6, lw = 1, linestyle = '-.', color=(0.3373, 0.6784, 0.8706, 1))
    plt.plot(d, (Ki3 + Ki2)/1e6, lw = 1, color='black', linestyle = '-')

    plt.axhline(y=0, color='black', linestyle='-', lw=0.4)
    plt.xlabel('Fracture depth ($m$)')
    plt.ylabel('Stress intensity factors ($MPa^{1/2}$)')
    plt.xlim(0, 200)
    plt.ylim(-20, 20)
    plt.tight_layout()

    if save_nice_fig:
      plt.savefig(f'K1-3_Rxx_{round(Rxx*1e-6,2)}.png', dpi = 400)

    plt.show()

else:
  #plot Fig. 10
  plt.plot(d/1000, Kic/1e6)
  plt.axhline(y=0.0, color='black', linestyle='--')
  #plt.axhline(y=0.1, color='blue', linestyle='--')
  plt.axhline(y=0.4, color='red', linestyle='-.')

  plt.xlabel('d')
  plt.ylabel('Kic')
  plt.xlim(0, 250)
  plt.ylim(-4, 12)
  #plt.gca().invert_yaxis()
  plt.show()




#unlooped version below
sys.exit()

#plot Fig. 4
plt.plot(d/H, flambda_out)
plt.xlabel('d/H')
plt.ylabel('F')
plt.xlim(0, 1)
plt.ylim(0, 25)

plt.show()
sys.exit()

#fix the function to values of d and H required
G_fixed = lambda b: G(b, d, H)

#calculate integral within bounds
integral, err = quad(G_fixed, 0, d)

#calculate stress intensity (Eq. 18)
Ri1 = 1.12*Rxx*( ( np.pi * d)**(1/2) )
Ri2 = ( ( 2*rhoi*g )/( ( np.pi * d)**(1/2) ) )*integral
Kic = Ri1 + Ri2

print(Kic/1e6)

